
library(countrycode)
library(foreign)
library(sandwich)
library(mice)
library(plyr)
library(VGAM)
library(nnet)
require(boot)
library(mlogit)
library(AER)
library(reshape2)
library(clusterSEs)
library(readstata13)
library(tidyr)
library(stargazer)
library(margins)
library(ggplot2)
library(extrafont)
options(warn=1)


###standard errors function to replicate stata (hainmueller code)
require(sandwich)
require(lmtest)
library(RColorBrewer)

library(plm)
library(MASS)
library(quantmod)
vcovCluster <- function(model, cluster)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)          
  K <- model$rank  
  if(M<50){
    warning("Fewer than 50 clusters, variances may be unreliable (could try block bootstrap instead).")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}


##########functions to automatically grab CIs etc from multinom

multinom_function <- function(data, choice1, choice2, iterations){
  nc_data<-data[complete.cases(data), ]
  ml_data<-mlogit.data(nc_data, shape='wide', choice=choice1, chid.var = NULL, alt.levels = NULL)
  resp <- names(ml_data)[1]
  covar <- setdiff(names(ml_data)[-1], c('country_num', 'chid', 'alt', choice2))
  ff <- as.formula(paste(resp, "~", "0","|", paste(covar, collapse="+")))
  Model <- mlogit(ff, ml_data)
  CIs<-cluster.bs.mlogit(Model,ml_data, ~as.factor(country_num), boot.reps=iterations)$ci
  
  est_weak<-Model$coefficients[3]
  lower_weak<-CIs[3, 1]
  upper_weak<-CIs[3, 2]
  weak_cis<-cbind(est_weak, lower_weak, upper_weak, '95%', 'base0', 'weak')
  colnames(weak_cis)<-c('est', 'lower', 'upper', "CI", 'base', 'type')
  
  est_strong<-Model$coefficients[4]
  lower_strong<-CIs[4, 1]
  upper_strong<-CIs[4, 2]
  strong_cis<-cbind(est_strong, lower_strong, upper_strong, '95%', 'base0', 'strong')
  colnames(strong_cis)<-c('est', 'lower', 'upper', "CI", 'base', 'type')
  
  
  Model90 <- mlogit(ff, ml_data)
  CIs<-cluster.bs.mlogit(Model90,ml_data, ~as.factor(country_num), ci.level=.9, boot.reps=iterations)$ci
  
  est_weak90<-Model90$coefficients[3]
  lower_weak90<-CIs[3, 1]
  upper_weak90<-CIs[3, 2]
  weak90_cis<-cbind(est_weak90, lower_weak90, upper_weak90, '90%', 'base0', 'weak')
  colnames(weak90_cis)<-c('est', 'lower', 'upper', "CI", 'base', 'type')
  
  est_strong90<-Model90$coefficients[4]
  lower_strong90<-CIs[4, 1]
  upper_strong90<-CIs[4, 2]
  strong90_cis<-cbind(est_strong90, lower_strong90, upper_strong90, '90%', 'base0', 'strong')
  colnames(strong90_cis)<-c('est', 'lower', 'upper', "CI", 'base', 'type')
  
  
  
  ####base1
  ml_data<-mlogit.data(nc_data, shape='wide', choice=choice2, chid.var = NULL, alt.levels = NULL)
  resp <- names(ml_data)[2]
  covar <- setdiff(names(ml_data)[-1][-1], c('country_num', 'chid', 'alt', choice1))
  ff <- as.formula(paste(resp, "~", "0","|", paste(covar, collapse="+")))
  Model <- mlogit(ff, ml_data)
  CIs<-cluster.bs.mlogit(Model,ml_data, ~as.factor(country_num), boot.reps=iterations)$ci
  
  est_strong_b1<-Model$coefficients[4]
  lower_strong_b1<-CIs[4, 1]
  upper_strong_b1<-CIs[4, 2]
  strong_b1_cis<-cbind(est_strong_b1, lower_strong_b1, upper_strong_b1, '95%', 'base1', 'strong')
  colnames(strong_b1_cis)<-c('est', 'lower', 'upper', "CI", 'base', 'type')
  
  
  Model90 <- mlogit(ff, ml_data)
  CIs<-cluster.bs.mlogit(Model90,ml_data, ~as.factor(country_num), ci.level=.9, boot.reps=iterations)$ci
  
  est_strong_b190<-Model90$coefficients[4]
  lower_strong_b190<-CIs[4, 1]
  upper_strong_b190<-CIs[4, 2]
  strong_b190_cis<-cbind(est_strong_b190, lower_strong_b190, upper_strong_b190, '90%', 'base1', 'strong')
  colnames(strong_b190_cis)<-c('est', 'lower', 'upper', "CI", 'base', 'type')
  
  df<-data.frame(rbind(weak_cis, strong_cis, weak90_cis, strong90_cis, strong_b1_cis, strong_b190_cis))
  df$est<-as.numeric(as.character(df$est))
  df$lower<-as.numeric(as.character(df$lower))
  df$upper<-as.numeric(as.character(df$upper))
  
  df
}

multinom_plot<-function(df){
  ggplot(df,
         aes(x=factor(paste(type, base), levels=c('strong base1', 'strong base0', 'weak base0')), y=as.numeric(as.character(est)), 
             group=factor(CI, levels=c('90%', '95%')), linetype=factor(CI, levels=c('90%', '95%'))))+
    geom_point()+
    theme_bw()+geom_pointrange(size=.5, ymin=as.numeric(as.character(df$lower)), 
                               ymax=as.numeric(as.character(df$upper)))+
    ylim((min(df$lower)-abs(.15*min(df$lower))), (max(df$upper)+.15*max(df$upper)))+scale_linetype_manual(values=c('solid', 'dotted'),
                                                                                                          guide = guide_legend(override.aes = list(
                                                                                                            linetype = c("solid", "dotted"),
                                                                                                            shape=c(NA, NA)) ))+
    theme(legend.position="bottom",  axis.title=element_text(size=16),
          text=element_text(family="Times", size=16),
          axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          legend.title=element_blank())+
    coord_flip() +
    geom_hline(aes(yintercept=0), lty=4) +
    ylab("") +
    xlab("")+scale_x_discrete(labels=c('Overt (base covert)', 'Overt (base none)', 'Covert (base none)'))
  
}

iterations=1000

##Conditional analysis DF [main analysis]
exprop_dem_gdp<-read.csv('exprop_dem_gdp_2021updated.csv', stringsAsFactors=F)

##Full country panel [main analysis]
ex_full_panel<-read.csv('ex_full_panel_2021updated.csv', stringsAsFactors=F)

##Panel of countries with FDI stock [appendix]
fdi_stock_panel<-read.csv('fdi_stock_panel_2021updated.csv', stringsAsFactors=F)

##ISDS [appendix]
isds<-read.csv('isds_merged.csv', stringsAsFactors=F)

##Only Kobrin countries [appendix]
ex_kobrin_panel<-read.csv('ex_kobrin_panel_2021updated.csv', stringsAsFactors=F)


################################################################################
####Paper Figure 1############################################################
################################################################################
##all expropriation events data

#defining overt expropriation
exprop_dem_gdp$strong_exprop<-NA
exprop_dem_gdp$strong_exprop[exprop_dem_gdp$type=='nat']<-1
exprop_dem_gdp$strong_exprop[exprop_dem_gdp$type=='cont']<-0
exprop_dem_gdp$strong_exprop[exprop_dem_gdp$type=='int']<-0
exprop_dem_gdp$strong_exprop[exprop_dem_gdp$type=='sale']<-0
###plotting expropriation events by year
strong_exprop_by_year<-ddply(exprop_dem_gdp, .(year), summarize, num_events=sum(strong_exprop, na.rm=T))
weak_exprop_by_year<-ddply(exprop_dem_gdp, .(year), summarize, num_events=length(strong_exprop[strong_exprop==0]))

strong_exprop_by_year$type<-'overt'
weak_exprop_by_year$type<-'covert'

exprop_by_year<-rbind(strong_exprop_by_year, weak_exprop_by_year)
exprop_by_year<-subset(exprop_by_year, is.na(exprop_by_year$year)==FALSE)

pdf('expropriation_events_by_year.pdf', width=8, height=6)
ggplot(exprop_by_year, aes(x=year, y=num_events, group=type, linetype=type))+
  geom_line()+theme_bw()+theme(axis.title=element_text(size=16),
                               text=element_text( size=16))+
  xlab('')+ylab('Expropriation Events per Year')+xlim(c(1960,2014))
dev.off()

#Number of overt vs. covert
nrow(subset(exprop_dem_gdp, strong_exprop==0))
nrow(subset(exprop_dem_gdp, strong_exprop==1))

###############################################################
#######Table 1 Conditional Legislative Constraints############
##############################################################
m1_lc<-glm(strong_exprop~leg_constraints, data=exprop_dem_gdp, family='binomial')
cov.m1_lc_cluster <- vcovCL(m1_lc, ~country)
robust.se.m1_lc_cluster <- sqrt(diag(cov.m1_lc_cluster))

m2_lc<-glm(strong_exprop~leg_constraints+democracy_p+ll_gdp_pc+ll_gdp_pc2+factor(sector2), data=exprop_dem_gdp, family='binomial')
cov.m2_lc_cluster <- vcovCL(m2_lc, ~country)
robust.se.m2_lc_cluster <- sqrt(diag(cov.m2_lc_cluster))

m3_lc<-glm(strong_exprop~leg_constraints+ll_gdp_pc+ll_gdp_pc2+factor(sector2)+region+factor(year_bin), data=exprop_dem_gdp, family='binomial')
cov.m3_lc_cluster <- vcovCL(m3_lc, ~country)
robust.se.m3_lc_cluster <- sqrt(diag(cov.m3_lc_cluster))

m4_lc<-glm(strong_exprop~leg_constraints+democracy_p+ll_gdp_pc+ll_gdp_pc2+factor(sector2)+region+factor(year_bin), data=exprop_dem_gdp, family='binomial')
cov.m4_lc_cluster <- vcovCL(m4_lc, ~country)
robust.se.m4_lc_cluster <- sqrt(diag(cov.m4_lc_cluster))

exp0<-(exp(m1_lc$coefficients[1]+quantile(exprop_dem_gdp$leg_constraints, .25, na.rm=T)*m1_lc$coefficients[2]))/(1+(exp(m1_lc$coefficients[1]+quantile(exprop_dem_gdp$leg_constraints, .25, na.rm=T)*m1_lc$coefficients[2])))
exp1<-(exp(m1_lc$coefficients[1]+quantile(exprop_dem_gdp$leg_constraints, .75, na.rm=T)*m1_lc$coefficients[2]))/(1+(exp(m1_lc$coefficients[1]+quantile(exprop_dem_gdp$leg_constraints, .75, na.rm=T)*m1_lc$coefficients[2])))
effects<-exp1-exp0

clustered_conditional_expropriation_legconstraints_table<-stargazer(m1_lc, m2_lc, m3_lc, m4_lc, se=list(robust.se.m1_lc_cluster, robust.se.m2_lc_cluster, robust.se.m3_lc_cluster, robust.se.m4_lc_cluster),
                                                                    align=TRUE, omit=c('region','year_bin','Constant'), header=FALSE, label="table:conditional_expropriation_legconstraints",
                                                                    covariate.labels=c("Legislative Constraints", "Democracy", "GDP per capita", "GDP per capita squared", "Extractive sector",
                                                                                       "Financial sector", "Manufacturing sector", "Services sector", "Utilities sector"), title="Legislative Constraints and Use of Overt Expropriation Conditional on Expropriating", font.size="footnotesize", column.sep.width="0pt",
                                                                    add.lines=list(c("Decade/Region FE", "N", "N", "Y","Y")),
                                                                    dep.var.labels="Use of overt (1) vs. covert (0) expropriation", no.space=TRUE,
                                                                    notes="Robust standard errors clustered at country level")
write(clustered_conditional_expropriation_legconstraints_table[4:length(clustered_conditional_expropriation_legconstraints_table)-1], 'R_analysis/PaperOutput/conditional_expropriation_legconstraints_table.tex')


###Interpretation of substantive effects
#finding most common region and sector
table(exprop_dem_gdp$sector2)
table(exprop_dem_gdp$region)
table(exprop_dem_gdp$democracy_p)
table(exprop_dem_gdp$year_bin)
quantile(exprop_dem_gdp$leg_constraints,c( .25,.75), na.rm=T)
quantile(ex_full_panel$leg_constraints,c( .25,.75), na.rm=T)

####substantive effects

prediction_lc<-with(exprop_dem_gdp, data.frame(leg_constraints=c(0,.476),
                                            democracy_p=0,ll_gdp_pc=mean(na.omit(exprop_dem_gdp$ll_gdp_pc)),
                 ll_gdp_pc2=mean(na.omit(exprop_dem_gdp$ll_gdp_pc2)), sector2="extractive",
                 region="afr",year_bin="1970s"))


predict(m4_lc, prediction_lc,type="response")




###############################################################
##############Table 2 Conditional Vertical Accountability########
##############################################################
m1_ac<-glm(strong_exprop~v2x_veracc, data=exprop_dem_gdp, family='binomial')
cov.m1_ac_cluster <- vcovCL(m1_ac, ~country)
robust.se.m1_ac_cluster <- sqrt(diag(cov.m1_ac_cluster))

m2_ac<-glm(strong_exprop~v2x_veracc+democracy_p+ll_gdp_pc+ll_gdp_pc2+factor(sector2), data=exprop_dem_gdp, family='binomial')
cov.m2_ac_cluster <- vcovCL(m2_ac, ~country)
robust.se.m2_ac_cluster <- sqrt(diag(cov.m2_ac_cluster))

m3_ac<-glm(strong_exprop~v2x_veracc+ll_gdp_pc+ll_gdp_pc2+factor(sector2)+region+factor(year_bin), data=exprop_dem_gdp, family='binomial')
cov.m3_ac_cluster <- vcovCL(m3_ac, ~country)
robust.se.m3_ac_cluster <- sqrt(diag(cov.m3_ac_cluster))

m4_ac<-glm(strong_exprop~v2x_veracc+democracy_p+ll_gdp_pc+ll_gdp_pc2+factor(sector2)+region+factor(year_bin), data=exprop_dem_gdp, family='binomial')
cov.m4_ac_cluster <- vcovCL(m4_ac, ~country)
robust.se.m4_ac_cluster <- sqrt(diag(cov.m4_ac_cluster))

conditional_expropriation_vertacc_table<-stargazer(m1_ac, m2_ac, m3_ac,m4_ac, se=list(robust.se.m1_ac_cluster, robust.se.m2_ac_cluster, robust.se.m3_ac_cluster,robust.se.m4_ac_cluster),
                                                   align=TRUE, omit=c('year_bin', 'region',"Constant"), header=FALSE, label="table:conditional_expropriation_vertacc",
                                                   covariate.labels=c("Vertical Accountability", "Democracy", "GDP per capita", "GDP per capita squared", "Extractive sector",
                                                                      "Financial sector", "Manufacturing sector", "Services sector", "Utilities sector"), title="Vertical Accountability and use of Overt Expropriation Conditional on Expropriating", font.size="footnotesize", column.sep.width="0pt",
                                                   add.lines=list(c("Decade/Region FE", "N", "N", "Y","Y")),
                                                  dep.var.labels = "Use of overt (1) vs. covert (0) expropriation", no.space=TRUE,
                                                   notes="Robust standard errors clustered at the country level")
write(conditional_expropriation_vertacc_table[4:length(conditional_expropriation_vertacc_table)-1], 'R_analysis/PaperOutput/conditional_expropriation_vertacc_table.tex')

####substantive effects

quantile(exprop_dem_gdp$v2x_veracc,c( .25,.75), na.rm=T)

prediction_ac<-with(exprop_dem_gdp, data.frame(v2x_veracc=c(-0.913,0.627),
                                               democracy_p=0,ll_gdp_pc=mean(na.omit(exprop_dem_gdp$ll_gdp_pc)),
                                               ll_gdp_pc2=mean(na.omit(exprop_dem_gdp$ll_gdp_pc2)), sector2="extractive",
                                               region="afr",year_bin="1970s"))


predict(m4_ac, prediction_ac,type="response")

###############################################################
#######Table 3 Panel Legislative Constraints############
##############################################################
#Creating decade bins
ex_full_panel$year_bin<-NA
ex_full_panel$year_bin[ex_full_panel$year>=1960&ex_full_panel$year<1970]<-"1960s"
ex_full_panel$year_bin[ex_full_panel$year>1969&ex_full_panel$year<1980]<-"1970s"
ex_full_panel$year_bin[ex_full_panel$year>1979&ex_full_panel$year<1990]<-"1980s"
ex_full_panel$year_bin[ex_full_panel$year>1989&ex_full_panel$year<2000]<-"1990s"
ex_full_panel$year_bin[ex_full_panel$year>1999&ex_full_panel$year<2010]<-"2000s"
ex_full_panel$year_bin[ex_full_panel$year>2009]<-"2010s"



m1_panel_lc_all_p<-lm(num_exprop~leg_constraints+democracy_p+log(lag_gdp_pc)+I(log(lag_gdp_pc)^2)+history+region+factor(year_bin), data=ex_full_panel)
m1_panel_lc_all_se_p<-sqrt(diag(vcovCluster(m1_panel_lc_all_p, ex_full_panel$country[is.na(ex_full_panel$democracy_p)==FALSE&is.na(ex_full_panel$leg_constraints)==FALSE&is.na(ex_full_panel$num_exprop)==FALSE])))

m1_panel_lc_strong_p<-lm(num_strong_exprop~leg_constraints+democracy_p+log(lag_gdp_pc)+I(log(lag_gdp_pc)^2)+history+region+factor(year_bin), data=ex_full_panel)
m1_panel_lc_strong_se_p<-sqrt(diag(vcovCluster(m1_panel_lc_strong_p, ex_full_panel$country[is.na(ex_full_panel$democracy_p)==FALSE&is.na(ex_full_panel$leg_constraints)==FALSE&is.na(ex_full_panel$num_strong_exprop)==FALSE])))

m1_panel_lc_weak_p<-lm(num_weak_exprop~leg_constraints+democracy_p+log(lag_gdp_pc)+I(log(lag_gdp_pc)^2)+history+region+factor(year_bin), data=ex_full_panel)
m1_panel_lc_weak_se_p<-sqrt(diag(vcovCluster(m1_panel_lc_weak_p, ex_full_panel$country[is.na(ex_full_panel$democracy_p)==FALSE&is.na(ex_full_panel$leg_constraints)==FALSE&is.na(ex_full_panel$num_weak_exprop)==FALSE])))

ols_panel_lc<-stargazer(m1_panel_lc_all_p, m1_panel_lc_strong_p, m1_panel_lc_weak_p, 
                        se=list(m1_panel_lc_all_se_p, m1_panel_lc_strong_se_p, m1_panel_lc_weak_se_p),
                        align=TRUE, omit=c('year_bin', 'region','Constant'), header=FALSE, label="table:ols_lc",
                        covariate.labels=c("Legislative Constraints", "Democracy", "GDP per capita", "GDP per capita squared", "Expropriation History"
                                           ), title="Legislative constraints and likelihood of using different types of expropriation, OLS panel", font.size="footnotesize", column.sep.width="15pt",
                        add.lines=list(c("Region/Year FE", "Y", "Y", "Y")),
                        omit.stat=c("f", "ser"),dep.var.caption  = "DV: Number of expropriations in given year",
                        dep.var.labels=c("Any expropriation","Overt expropriation","Covert expropriation"), no.space=TRUE,
                        notes = "Robust standard errors clustered at the country level")

write(ols_panel_lc, 'R_analysis/PaperOutput/ols_panel_lc.tex')

###############################################################
#######Table 4 Panel Vertical Accountability############
##############################################################

m1_panel_ac_all_p<-lm(num_exprop~v2x_veracc+democracy_p+log(lag_gdp_pc)+I(log(lag_gdp_pc)^2)+history+region+year_bin, data=ex_full_panel)
m1_panel_ac_all_se_p<-sqrt(diag(vcovCluster(m1_panel_ac_all_p, ex_full_panel$country[is.na(ex_full_panel$democracy_p)==FALSE&is.na(ex_full_panel$v2x_veracc)==FALSE&is.na(ex_full_panel$num_exprop)==FALSE])))

m1_panel_ac_strong_p<-lm(num_strong_exprop~v2x_veracc+democracy_p+log(lag_gdp_pc)+I(log(lag_gdp_pc)^2)+history+region+year_bin, data=ex_full_panel)
m1_panel_ac_strong_se_p<-sqrt(diag(vcovCluster(m1_panel_ac_strong_p, ex_full_panel$country[is.na(ex_full_panel$democracy_p)==FALSE&is.na(ex_full_panel$v2x_veracc)==FALSE&is.na(ex_full_panel$num_strong_exprop)==FALSE])))

m1_panel_ac_weak_p<-lm(num_weak_exprop~v2x_veracc+democracy_p+log(lag_gdp_pc)+I(log(lag_gdp_pc)^2)+history+region+year_bin, data=ex_full_panel)
m1_panel_ac_weak_se_p<-sqrt(diag(vcovCluster(m1_panel_ac_weak_p, ex_full_panel$country[is.na(ex_full_panel$democracy_p)==FALSE&is.na(ex_full_panel$v2x_veracc)==FALSE&is.na(ex_full_panel$num_weak_exprop)==FALSE])))

ols_panel_ac<-stargazer(m1_panel_ac_all_p, m1_panel_ac_strong_p, m1_panel_ac_weak_p, 
                        se=list(m1_panel_ac_all_se_p, m1_panel_ac_strong_se_p, m1_panel_ac_weak_se_p),
                        align=TRUE, omit=c('year_bin', 'region',"Constant"), header=FALSE, label="table:ols_ac",
                        covariate.labels=c("Vertical Accountability", "Democracy", "GDP per capita", "GDP per capita squared", "Expropriation History"), 
                        title="Vertical Accountability and likelihood of using different types of expropriation, OLS panel", font.size="footnotesize", column.sep.width="15pt",
                        add.lines=list(c("Region/Decade FE", "Y", "Y", "Y")),dep.var.caption  = "DV: Number of expropriations in given year",
                        dep.var.labels=c("Any expropriation","Overt expropriation","Covert expropriation"),
                        omit.stat=c("f", "ser"),
                        no.space=TRUE,
                        notes = "Robust standard errors clustered at the country level")

write(ols_panel_ac, 'R_analysis/PaperOutput/ols_panel_ac.tex')


#######################################################################
#########Figure 2: multinomial logit - Legislative Constraints#########
########################################################################
#no controls
nc_lc<-multinom_function(ex_full_panel[c('exprop_numeric', 'exprop_numeric_b1', 'leg_constraints', 'country_num')], 'exprop_numeric', 'exprop_numeric_b1', iterations)
write.csv(nc_lc, 'nc_lc.csv')
ggsave('R_analysis/PaperOutput/multinom_nocontrols_lc.pdf', plot=multinom_plot(nc_lc), height=5, width=5.25)


#full controls + decade fes

fc_lc_df_decade_fe<-multinom_function(ex_full_panel[c('exprop_numeric', 'exprop_numeric_b1','leg_constraints','year_bin', 'democracy_p', 
                                              'll_gdp_pc', 'll_gdp_pc2', 'history', 'region', 'country_num')], 'exprop_numeric', 'exprop_numeric_b1', iterations)

write.csv(fc_lc_df_decade_fe, 'fc_lc_df_decade_fe.csv')
ggsave('R_analysis/PaperOutput/multinom_fullcontrols_lc_plot_decade_fe.pdf', plot=multinom_plot(fc_lc_df_decade_fe), height=5, width=5.25)



#full controls - decade fes
#fc_lc_df<-multinom_function(ex_full_panel[c('exprop_numeric', 'exprop_numeric_b1','leg_constraints', 'democracy_p', 
#                                                      'll_gdp_pc', 'll_gdp_pc2', 'history', 'region', 'country_num')], 'exprop_numeric', 'exprop_numeric_b1', iterations)
#write.csv(fc_lc_df, 'fc_lc_df.csv')
#ggsave('R_analysis/PaperOutput/multinom_fullcontrols_lc_plot.pdf', plot=multinom_plot(fc_lc_df), height=5, width=5.25)


#######################################################################
#########Figure 3: multinomial logit - Vertical Accountability#########
########################################################################
#no controls
nc_ac_df<-multinom_function(ex_full_panel[c('exprop_numeric', 'exprop_numeric_b1', 'v2x_veracc', 'country_num')], 'exprop_numeric', 'exprop_numeric_b1', iterations)
write.csv(nc_ac_df, 'nc_ac_df.csv')
ggsave('R_analysis/PaperOutput/multinom_nocontrols_ac_plot.pdf', plot=multinom_plot(nc_ac_df), height=5, width=5.25)

#full controls + decade fes
fc_ac_df_decade_fes<-multinom_function(ex_full_panel[c('exprop_numeric', 'exprop_numeric_b1', 'v2x_veracc','year_bin', 'democracy_p',
                                            'll_gdp_pc', 'll_gdp_pc2', 'history', 'region', 'country_num')], 'exprop_numeric', 'exprop_numeric_b1', iterations)

data<-ex_full_panel[c('exprop_numeric', 'exprop_numeric_b1', 'v2x_veracc','year_bin', 'democracy_p',
                      'll_gdp_pc', 'll_gdp_pc2', 'history', 'region', 'country_num')]

nc_data<-data[complete.cases(data), ]
choice1='exprop_numeric'
choice2='exprop_numeric_b1'
ml_data<-mlogit.data(nc_data, shape='wide', choice=choice1, chid.var = NULL, alt.levels = NULL)
resp <- names(ml_data)[1]
covar <- setdiff(names(ml_data)[-1], c('country_num', 'chid', 'alt', choice2))
ff <- as.formula(paste(resp, "~", "0","|", paste(covar, collapse="+")))
Model <- mlogit(ff, ml_data)
CIs<-cluster.bs.mlogit(Model,ml_data, ~as.factor(country_num), boot.reps=iterations)$ci

####base1
ml_data<-mlogit.data(nc_data, shape='wide', choice=choice2, chid.var = NULL, alt.levels = NULL)
resp <- names(ml_data)[2]
covar <- setdiff(names(ml_data)[-1][-1], c('country_num', 'chid', 'alt', choice1))
ff <- as.formula(paste(resp, "~", "0","|", paste(covar, collapse="+")))
Model <- mlogit(ff, ml_data)
CIs<-cluster.bs.mlogit(Model,ml_data, ~as.factor(country_num), boot.reps=iterations)$ci



write.csv(fc_ac_df_decade_fes, 'fc_ac_df_decade_fes.csv')
ggsave('R_analysis/PaperOutput/multinom_fullcontrols_ac_plot_decade_fes.pdf', plot=multinom_plot(fc_ac_df_decade_fes), height=5, width=5.25)


#full controls
#fc_ac_df<-multinom_function(ex_full_panel[c('exprop_numeric', 'exprop_numeric_b1', 'v2x_veracc', 'democracy_p',
#                                              'll_gdp_pc', 'll_gdp_pc2', 'history', 'region', 'country_num')], 'exprop_numeric', 'exprop_numeric_b1', iterations)
#write.csv(fc_ac_df, 'fc_ac_df.csv')
#ggsave('R_analysis/PaperOutput/multinom_fullcontrols_ac_plot.pdf', plot=multinom_plot(fc_ac_df), height=5, width=5.25)

